function [elements,processInfo] = ...
    CalculateMaterialResponse(...
            elements,model,processInfo)
ne = model.countElements();
dim = 2;
for i=1:ne
    e = model.getElement(i-1);
    f = e.getFace();
    ue = f.getUHatLocal();
    nn = f.countModes();
    ng = e.getNumGaussPoint();
    Ce = elements(i).elasticConstitutiveMatrix;
    B = zeros(3,dim*nn);
    for j=1:ng
        xi = e.getGaussPoint(j-1);
        %compute strain
        B = e.computeB(xi(1),xi(2),f,B);
        epsilon = B*ue;
        sigmae = Ce*epsilon;
        if strcmp(processInfo.integrationType,'implicit')
            %compute equivalent strain
            etilde = computeEquivalentStrain...
                (elements(i),epsilon);
            %check damage function and update internal variable
            if etilde > elements(i).kappa_old(j)
                %loading
                elements(i).kappa(j) = etilde;
                isLoad = 1;
            else
                %unloading
                elements(i).kappa(j) = ...
                    elements(i).kappa_old(j);
                isLoad = 0;
            end
            %compute damage variable
            d = computeDamageVariable(elements(i),...
                elements(i).kappa(j));
            %compute stress
            elements(i).stress(:,j) = (1-d)*sigmae;
            %compute algorithmic tangent
            if isLoad
                alpha = elements(i).kappa(j); %to make consistent with update scheme by Oliver
                tmp = computeDamageSlope...
                    (elements(i),alpha)/alpha;
                if strcmp(elements(i).EquivalentStrainMode,...
                        'euclidean')
                    epsilon_noneng = [epsilon(1);...
                        epsilon(2);epsilon(3)/2];
                    elements(i).constitutiveMatrix(:,:,j) =...
                        (1-d)*Ce - ...
                        tmp*(epsilon_noneng*sigmae');
                elseif strcmp(elements(i).EquivalentStrainMode,...
                        'strain norm')
                    elements(i).constitutiveMatrix(:,:,j) =...
                        (1-d)*Ce - ...
                        tmp*(sigmae*sigmae');
                end
                %             alpha
                %             d
                %             epsilon_noneng
                %elements(i).constitutiveMatrix(:,:,j)
            else
                elements(i).constitutiveMatrix(:,:,j) =...
                    (1-d)*Ce;
            end
            elements(i).damage(j) = d;
            if(abs(d-1)<1e-6)
                elements(i).kappa(j)
                i
            end
        elseif strcmp(processInfo.integrationType,'implex')
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % perform an implicit update
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %compute equivalent strain
            etilde = computeEquivalentStrain...
                (elements(i),epsilon);
            %check damage function and update internal variable
            if etilde > elements(i).kappa_old(j)
                %loading
                elements(i).kappa(j) = etilde;
            else
                %unloading
                elements(i).kappa(j) = ...
                    elements(i).kappa_old(j);
            end
            %compute damage variable
            d = computeDamageVariable(elements(i),...
                elements(i).kappa(j));
            elements(i).damage(j) = d;
            elements(i).stressImplicit(:,j) = ...
                (1-d)*sigmae;
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % perform an explicit update
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %compute extrapolated strain-like variable
            kappan = elements(i).kappa_old(j);
            kappanm1 = elements(i).kappa_old_old(j);
            if processInfo.Dt_old == 0
                tmp = 0;
            else
                tmp = processInfo.Dt/processInfo.Dt_old;
            end
            if tmp < 0
                tmp = 0;
            end
            kappa = kappan + tmp*(kappan - kappanm1);
            %compute damage variable based on extrapolated strain-like
            %variable
            d = computeDamageVariable(elements(i),...
                kappa);
            %compute extrapolated stress
            elements(i).stress(:,j) = (1-d)*sigmae;
            %compute extrapolated algorithmic tangent
            elements(i).constitutiveMatrix(:,:,j) =...
                        (1-d)*Ce;
        end
    end
end
end
